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ABSTRACT 



The development of microelectronics has brought into being 
large capacity digital memories in a small package. In the 
foreseeable future even more advances can be seen in this 
trend. Therefore the use of digital computers in control sys- 
tems will play an even larger role than today. 

This work involves a fourth order system to simulate the 
control and dynamics of a missile. Proportional navigation 
is used as the guidance method. Studied are the effects of 
applying different controls which are considered best from a 
computer study and the effects of applying digital filtering 
methods. Although these studies were applied to a specific 
problem, an attempt is made to keep the discussion general 
in order that the methods may be considered for other problems . 
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CHAPTER I 
INTRODUCTION 

If it is desired to hit a target with a missile, some form of 
guidance must be employed. This guidance would be selected 
weighing the desired objectives and an estimate of the target 
capabilities against the amount of money and missile space 
available. As suggested in the introduction, the advances in 
computer size and capabilities is, and will in the future, con- 
tinue to change the weighting of these factors . Thus in the 
foreseeable future it may be possible to design a model of the 
system which will be in sufficient detail to predict the respon- 
se of the real system to a high degree of accuracy. Using the 
micro-computers, this model could then become a part of the 
control system in the missile. Such a control system might be 
described in the block diagram in Figure 1-1 . 

In Figure 1-1 the missile measures and reacts to a signal 
which depends on the type of guidance employed. Such a 
signal might depend on the rate of change in the missile's 
radar seeker head which is tracking the line of sight from the 
missile to target. The missile will react to this signal, pro- 
ducing through its dynamics a change in the direction of its 
velocity vector, or the missile will take some other action 
called for by the particular guidance law being used. The 
model receives the measured signals and considers the 
measurable states from the missile dynamics. Since all the 
states of the model are measurable and the digital model in 
itself generates relatively little noise, the model becomes a 
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Figure 1-1. Block Diagram of a Dynamic System. 
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state predictor or estimator. 

Since all the states of the model are available, a control 
law considering all the states can now be designed. This 
control law is applied in the controller to guide the missile. 

If the desired objectives of the system can justify the 
above approach, the system performance of the modified system 
should be an improvement over most conventional guidance 
systems . 
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CHAPTER II 

PROPORTIONAL NAVIGATION 

> 

Adler in [9)] defined proportional navigation as, 

A course in which the rate of change of missile heading is 
directly proportional to the rate of rotation of the line-of- 
sight from the missile to the target. 

y= Kg (1) 

where 

y is the angular rate of change of the missile velocity vector, 
a is the angular rate of change of the line-of- sight. 

K is a constant, typically between three and five. 

To gain a better understanding of proportional navigation, 
some of the simpler forms of line-of-sight guidance will be 
examined. Pursuit course (sometimes called pure pursuit) always 
aims the missile directly at the target along the line-of-sight. 
This method will always end in a tail chase, even if the missile 
is launched head-on with the target. High missile acceler- 
ations are required. It can be shown from the equations of 
motion that if the missile to target speed ratio^ exceeds two, 
the final missile turning rate will approach infinity. Thus 
pursuit course guidance may not be very satisfactory although 
it is very simple to implement. 

The next step is to lead the line-of-sight by an angle 
which is a function of the target velocity. This type is called 
constant bearing course, and is achieved by aligning the 
relative missile to target velocity vector with the line-of- 
sight. That is to say the line-of-sight maintains a constant 
direction in space. It can be shown that this method can 
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cause a missile to follow a target to a collision even if the 
target is maneuvering. This method however requires an in- 
stantaneous correction to each change in the line-of-sight. 

Consider the constant bearing trajectory shown in 
Figure 2-1. In the figure the line-of-sight is shown for 
several samples made by the radar. As long as the missile 
remains on the collision course, the angular rate of change of 
the line-of-sight will remain zero and no control action is 
required. Any rotation of the line-of-sight, indicating a de- 
parture from the collision course will be detected by the 
missile's radar. The missile using proportional pavigation 
will turn at a rate proportional to this rotation in a direction 
to reduce the line-of-sight rate and return to a constant 
bearing collision course. 

In the Figure 2-2 above, the lead angle /3is defined: 
p=o-y (2) 

This angle (/3) is independent of any reference used to define 
ct and y. If, as in the discussion above, we could achieve a 
constant bearing course, £ would remain a constant, and 
therefore ^ would remain zero. 

Thus for a constant bearing course , 

($= a - y = 0 (3) 

Since the missile cannot react instantaneously, /}is not 
always zero. Thus for any instant of time (1) may be written 
as , 

K,a= 0 + y (4) 

9 

If /3 is equal to zero, the missile is on a constant bearing 
collision course . Therefore the object of the control system 
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Figure 2-1. Constant Bearing Trajectory. 
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Figure 2-2. Proportional Navigation Geometry. 
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Figure 2-3. Missile Guidance Inter Loop. 
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Figure 2-4. Missile Guidance Outer Loop. 
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in proportional navigation is to minimize fiand if a becomes a 
constant or zero to drive |8to zero. 

Guidance and Control 

A typical homing missile inter guidance loop is shown in 
Figure 2-3. The following discussion will be limited to only 
coplanar motion of the missile and target. A constant missile 
and target speed is also assumed throughout the rest of this 
paper. 

In Figure 2-3, the radar tracks the target using a servo 
loop to keep the antenna on the line-of-sight. A delayed sig- 
nal which is a measure of a, is obtained from the radar system. 
A voltage which is proportional to the rate of change of this 
signal is obtained. This voltage is compared with the output 
of the dynamic system y, and the error signal is obtained. 

This error signal is sent either to the hydraulic system to 
rotate the control surfaces or to a hydraulic/valve system 
which controls side thrust jets, depending on the type of 
missile being used. 

Kinematics 

As might be suggested by the name inter loop, there is 
also an outer loop as shown in Figure 2-4. In the Figure, 
the guidance and control block is shown as in Figure 2-3. 

The effect of yon the line-of-sight and the effect of the tar- 
get motion on the line-of-sight are represented in the missile- 
target kinematics block. 

This outer loop is not as accessible to measurement 
and control, but the understanding of its effects is necessary 
in order to perform a simulation of the entire system . From 
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the geometry of Figure 2-2, the following equations may be 
written: 

R = V t COS(o) - V m COS(/?) (5) 

R a = - V t SIN ( a) + V m SIN ( /3) (6) 

where 

V t = Target speed (assumed Constant) 

V = Missiie speed (assumed Constant) 
m 

a - y 

R = the line-of-sight distance 
R = the rate of change of the line-of-sight. 

Then equations (5) and (6) may be rewritten: 

R = V t COS(a) - V m COS( a - y) (7) 

Ra = -V t SIN (a) + V m SIN( a - y) (8) 

Differentiating (8) with respect to time: 

Ra+Rc= -V t COS( CT ) a+ V m COS( CT - y)(a- y) 

= -Ro- v m cos (o-y)y (9) 

Therefore 

2Ra+ Ra= -V m COS(a- y) y (10) 

let 



R = -R/T 



and 

V r =tR| 

where T is the time to go to the target. Then equation (10) 
may be represented by: 

(-2 R/T + RS) a = -V m COS(a - y) y (11) 



thus 

o/Y= -V m COS(o-y) (12) 

R/T (TS - 2) 

Where V m K COS(a- y) is called the effective navigation constant. 
R/T 
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Since V m is constant and variations in COS(or- y) and V r are 
small, this term is assumed constant and values of four to 
six are considered best. 
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CHAPTER III 
FILTER AND CONTROL 

Plant Description 

A linear, time -invariant dynamic system is described in 
the flow graph, Figure 3-1. 

The transfer function of Figure 3-1 may be stated: 

x out ( s > = V 0 " 1 + ; ♦ . ♦ ♦ + b 2 s + b t 
Xj n (s) S n + 3 n s 11 ^ + . . + 32 s + 

Equation (1) may be written in one differential equation using 
matrix notation . 

X = F X + D U, C = B X (2) 

Where 

X is the state vector (n x 1) 

U is a scalar of the system inputs and controls 
F is a matrix of constants. For the system described in 
Figure 3-1, F is defined in equation (3). 
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D is a matrix of constants (n x 1) described in equation 
(4) for the system in Figure 3-1. 
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Figure 3-1. Flow Diagram of the system. 
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0 

0 



( 4 ) 



D = 

0 

bJ 

B is an (1 x n) row vector, described in equation (5) for 
the system described in Figure 3-1 . 



B = [b 1 b 2 . . . b n ] (5) 

The solution to (2) is given by: 

X(t) = $(t - t Q ) X(t Q ) + (t - t x ) D dt x U(t Q ) 1 (6) 



Where 

4>(t - t Q ) = -C -1 (SI - F) 



Solving the differential equation (2) for the discrete solution. 



X(k + 1) = $>X(k) + DEL U(k) 


( 7 ) 


Where 




$= e^fr *o ^ 


( 8 ) 


DEL = J' t to e F ^ D dt x 


( 9 ) 


Note that k(t - t Q ) is represented by k. 





Titus, in reference [10] developed a digital computer 
program called PHIDEL. This program provides a solution to 
equations (8) and (9) and will be used as a subroutine to 
obtain $and DEL. A listing of the PHIDEL program is provided 
in Appendix one . 

Plant Control 

It is desired to design a set of feedback coefficients to 
control the plant in accordance with a selected performance 



1 U(t) is held constant over the interval (t - t Q ) and is 
equal to U(t Q ) . 
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index. This set of coefficients will modify the plant states and 
in effect will cause a movement of the Z-plane poles. This 
effect must be weighed against the feedback gains derived for 
the desired performance index. 

Let J(N) be defined as the performance index of a discrete 
sample-data system which will be minimized with respect to 
U(K), the control for the plant. 

N + , 

J(N) = Min E [X 1 - (K) Q X(K) + r JJ Z (K - 1)] (10) 

U(K) K = M 

Q is a (n x n) positive-definite symmetric constant matrix, 
and R is a positive scalar constant. 

The principle of optimality states that any portion of an 
optimal trajectory is also an optimal trajectory. Therefore 
equation (1) may be rewritten: 

J(N) = Min [ X t (N) Q X(N) + r U 2 (N - 1) + J(N - 1)] (11) 

U(N) 

If the gradient of J(N) is taken with respect to U(N - 1) and 
set equal to zero, then an optimal control, U°(N - 1) may be 
determined. It is noted that J(N - 1) is not a function of 
U(N - 1). 

Using the above arguments, Titus developed the following 
algorithms: 

At (k + 1) = -jt P(k) 4- 

Atp(k) A+ r (12) 

P(k) = (k)P(k - 1) TO + Q + rA(k)A t (k) (13$ 

TO = ^+AAt(k) (14) 

Using the recursive relations of (12), (13) and (14), Titus [10] 
developed program OPTCON. This method is also discussed 
in references [3] and [4] by Titus, Strum and Demetry, and 
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in reference [6 ]by Ogden. A fortran listing of the program may 
be found in Appendix one . 

Digital Filter 

In the control discussion, the assumption was made that 
all the states of the dynamic system are observable states. This 
condition is certainly not always true and often the observable 
states may be measured only at the cost of some ambiguity due 
to the measurement noise. The inputs to the system such as 
the radar seeker discussed in Chapter Two will often be re- 
ceived in a noisy environment. 

Kalman in references [l] and [2] discussed these prob- 
lems and presented the theory for the desired filter. Schmidt 
03] ; Titus, Demetry and Strum [4] and Jardine [7] have dis- 
cussed this problem and developed methods of implementing this 
problem on the digital computer. In an effort to maintain clafity, 
some of these developments will be presented. 

The digital filter will provide a best estimate of all the 
states by weighing the past information with the present obser- 
vable states. This weighting is performed with the knowledge 
of the environment (excitation) and the measurement noise . 
Although in many cases the noise is very difficult to describe, 
the assumption that white noise is present is in general true. 
Thus if white noise can be discriminated, the integrity of the 
measured signals will be improved. 

Filter Design 

For a single variable, the desnity function f(x) for a 
gaussian distribution is given by: 

f (x) = 1 e'& x - u > 2 

l2ff)l a (15) 
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where 



£ 9 (x)f (x)dx = E [ g (x) ] (16) 

where E implies expected value, 
and 

^~f(x)dx = 1 (17) 

Assume that the systems measurement and excitation 
noise are white noise, where white noise has a gaussian dis- 
tribution and is spread uniformly over all frequencies. 

Noting that the input, U(k) to the plant is independent of 
U(k - 1), X(k) is independent of X(k - 1) and using the joint 
probability property of random independent samples, a density 
function for the state system may be written: 

F(X) = 1 e -£[X - X P" 1 [X - X ] ( 18 ) 

(2 vW 2 pi 

where X(k) is the best estimate of the states based on past 
information . 

Z(k - 1) = H X(k - 1) + V(k - 1) (19) 

Z is the noisy observable matrix, and V is the measurement 
noise vector defined by: 

E [V ] = 0 (20) 

E R (21) 

R is the covariance matrix of the measurement noise and P is 

the covariance matrix of the error. P is a symmetric matrix, 

(P„ = ) . The diagonal terms of P are equal to the variance 

of each state ( a ^ ) and the off-diagonal terms are equal to 

x i 

correlation between the states, 

(E[x t • x j ] - x . .£.). 

P = E[ [X - X]. [X - X] 1 ] 
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( 22 ) 



The expected value of the error is defined as the loss function: 

L = J (X - X) l (X - X)F(X/Z,X)dX (23) 

/ i A 

Taking the gradiant of equation (23) with respect to X; 

Vx L = J 2XF(X/Z,X)dX - 2X jF(X/Z,X)dX (24) 

Applying equations (16) and (17) and setting the result equal 
to zero 

V£L = 2 E X/Z,X - 2 X = 0 
Thus: 

X* = Max$) = E [ X/Z,X ] (25) 

X* is the best estimate of the states given the present noisy 
observable and the past prediction of the states. 

Filter Equations 

There are two methods of finding the recursive relations for 
the filter program. Method one assumes the random varables are 
gaussian and makes use of Bayes equation to derive an ex- 
pression for equation (25). Method two assumes a linear fil- 
ter, selects an algorithm for equation (25) and proves that this 
selection provides the best values for X*. 

Method One 

The covariance matrix of the observable error may be written: 

P z = E[(Z - Z)(Z - Z)*] (26) 

Combining equations (19) and (26), 

P z = E [ (HX + V - HX) (HX + V - HX)*] 

= HE[(X-X)(X-X) l ] + E [ W* ] 

Using equations (20), (21) and (22), 

P 2 = HPHt + R (27) 
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In order to find an expression for X*, Bayes formula is 

applied to equation (25) . 

X* = F(X/Z,X) = F(Z.X/X) F(X) 

F(Z,$) 

A 

Since Z(k) is independent of X(k) , equation (28) may be 
written: 

F (X/Z , X) = F (Z/X) F(X/X) F(X) 

F (Z) F(X) 

A a A 

Since F(X/X) = F(X/X) F(X) , equation (29) becomes, 

F (X) 

A A 

F(X/Z,X) = F (Z/X) F (X/X) 

F(Z) 

A 

Since X(k) is independent of X(k) , (30) becomes, 

F (X/Z,X) - F (Z/X) F (X) 

F(Z) 

Following the form in equation (18), F(Z) and F(Z/X) are 
defined: 

e - 



F(Z) = 1 e"^ z ” Z ) tr>z 1 ^ z z ) 

_-^vt R' 1 V) 



F (Z/X) = F(V) = 



1 



(2(ff.) n / 2 |Rl 2 
Combining equations (18), (31) and (33), 

A _1 p 

F(X/Z,X) = A e 
Where 



A = (HPH t + Rjg 

(277) n / 2 |r| i |P| i 



and 



( 28 ) 



(29) 



(30) 



(31) 



(32) 

(33) 



B = - (X - X) t P _1 (X - X) + (Z - Z) t (HPH t + R) -1 (Z - Z) 

ng the gradiant of B with respect to X equal to zero. 

V X B = -^(X^X^P _1 - 2(Z-Z) t (HPH t +R)" 1 '7 x (Z-Z) = 0 



since 



V (Z-Z) =^7y(HX+ V) = H, the algorithm, equation (35) 
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( 35 ) 



may be written. 

X* = X + PH^HPH 1 + R) -1 (Z - Z) 

Let the filter gain matrix G be defined . 

G = PH 1 (HPH 1 + R) -1 (36) 

Then the algorithm, equation (35) may be written: 

X* = X + G (Z - Z) (37) 

Method Two 

An alternate method of deriving the recursive relation for 
the gain matrix assumes a linear filter. The trace of the co- 
variance matrix is given by: 

L = E [ (X - X) t (X - X)] (38) 

It will be shown that if L is minimized with respect to the gain, 
then the gain will be given by equation (36). 

Substituting equation (37) into equation (38) , 

L = E [ (X-X) (X-X) 1 ] - E [ (X-X) (Z-Z^G 1 ] - 

E[G(Z-Z)(X-X) t ] + E [ G (Z-Z) (Z-Z) t G t ] (39) 

Since Z = HX + V 
Then Z = E [ HX + V ] = HX 

Assuming that E [(X-XjV^jis equal to zero, equation (39) may 
be written: 

L = E [ (X-X) (X-X) 1 ] - E [ (X-X) (X-X) 1 ] H^G 1 + GE [ W 1 ] 

- GH E [ (X-X) (X-X) 1 ] + GHE [ (X-X) (X-J^H k* 1 (40) 

Substituting equations (21) and (22) into equation (40), 

L = P - PH^ 1 - GHP + GHPH 1 ^ 1 + GRG 1 (41) 

Taking the gradiant of L with respect to G and solving for G. 

G = PH t (HPH 1 + R) -1 
Filter Programs 

Jardine, Titus, Demetry and Strum have designed programs 
to calculate the filter gains and the covariance matrix of error. 

29 



Such a program, (Program Filter) is listed in Appendix one. 
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CHAPTER IV 
DIGITAL SIMULATION 

If a digital filter can be obtained, then a method must be 
developed to evaluate this system. This chapter will describe 
the simulation of the missile tracking system, the digital pre- 
dictor, the control, and the kinematics. The results of this 
simulation for a particular missile will be discussed and tra- 
jectories from computer runs will be presented in Chapter Five. 
Missile Guidance 

Figure 4-1 is a block diagram of the digital simulation. As 
described in previous chapters, the missile system tracks the 
target, obtaining an input for the missile guidance. This signal 
is proportional to the angular rate of change of the line-of-sight 
angle. In this simulation it is assumed that white noise is 
added to this signal. 

The details of the missile guidance simulation are des- 
cribed in Appendix Two. This portion of the simulation de- 
scribes the tracking, control and dynamics of the missile as 
a transfer function relating the line-of-sight angular rate to the 
missile flight path angular rate. This transfer function is then 
formed into an F and D matrix as described in Chapter Three. 
Program Phidel is employed to obtain the <fand A matrixes for 
the state difference equation. 

The output of the missile guidance (y) is considered the 
observable for the predictor. Gaussian measurement noise 

with standard deviation (a v ) is added to the observable. The 

* 

output (y) is also an input to the kinematics. Figure 4-2 
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Figure 4-1 ♦ Block Diagram of Missile Simulation. 
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Figure 4-2. Discrete Flow Graph of the Plant. 
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describes the digital simulation of the guidance. The following 
fortran listing updates the discrete difference equation as dem- 
onstrated in Figure 4-2 . 

C THIS SECTION UPDATES THE STATE VECTOR X 
CALL RNDEV(NUNIF,DEV) 

W=SIGW*DEV 

CALL PROD (PHI ,X,KN ,KN , 1 , TEMPI) 

DO 803 I = 1 ,KN 
803 TEMP2 (1,1) = W*DEL(1 , 1 ) 

CALL ADD(XS,DINP,KN, 1 ,TELP) 

CALL PROD(AT ,TELP , 1 # KN , 1 ,TELP1) 

CALL PROD(DEL,TELPl ,KN, 1 , 1 ,TELP2) 

CALL ADD (TEMPI ,TEMP2 ,KN , 1 ,X) 

CALL ADD(X,TELP2 ,KN, 1 ,X) 

Digital filter-predictor 

The digital filter is similar to the filter described in 
reference [4 ]. The gain matrix (G) is evaluated each sample 
instead of using the steady-state values of the gain. The dis- 
crete flow graph of the filter is shown in Figure 4-3. 

One input to the filter is the noisy observable. This is 
the sum of the angular rate of the missile flight path (y) and 
the measurement noise. This noise is assumed gaussian with 
mean zero and variance a v . This measurement on the missile 
could be made by use of rate gyros. 

The other input to the filter is the deterministic input. 

* 

Here this would be the angular rate of the line-of-sight (ct) 
and the best estimate of its derivatives. In the simulation 
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Figure 4-3. Discrete Flow Graph of the Filter/predictor 
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these are obtained from the kinematics. In the missile this 
signal would be measured by the tracking system. Note that the 
present value of the deterministic input DI(k) is used as the 
input to the plant while the past value of the input DI(k - 1) is 
used in the filter. The reasoning is that the filter is taking the 
present value of the observable Z(k) and therefore must use the 
past value of the deterministic input DI(k - 1) which caused the 
observable Z(k). 

The fortran statements for the filter are listed below. 

C THIS SECTION CALCULATES XS , THE BEST ESTIMATE 
C OF THE STATE VECTOR 

CALL PROD(H,X # KP ,KN , 1 ,Y) 

DO 10 1=1 ,KP 

CALL RNDEV(NUNIF ,DEV) 

10 V(I,1) = SIGV(I)*DEV 
CALL ADD(Y,V,KP,1 ,Z) 

CALL PROD (PHI ,XS,KN,KN,1 , TEMPI) 

DO 11 I = 1 ,KN 
DO 11 J = 1,KN 

11 TEMP2 (I , J) = -TEMP2 (I , J) 

CALL PROD(TEMP2 ,TEMP1 ,KN,KN,1 ,TEMP3) 

CALL ADD(TEMP1 ,TEMP3 ,KN , 1 ,TEMP3) 

CALL ADD(XS ,DINP ,KN , 1 ,TELP) 

CALL PROD(AT,TELP, 1 ,KN,1 ,TEMP1) 

CALL PROD(D EL, TEMPI ,KN,1 ,1 ,TELP) 

CALL PROD(TEMP2 ,TELP,KN,KN, 1 ,TELP1) 

CALL ADD(TELP ,TELP1 ,KN , 1 ,TELP1) 

CALL PROD(G ,Z,KN ,KP, 1 ,TELP2) 
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CALL ADD (TEM P3 , TELP1 , KN , 1 , XS) 

CALL ADD(XS , TELP2 , KN , 1 , XS) 

The Kinematics 

The kinematics combines the flight path angular rate of the 
missile, the missile speed (assumed constant) and the target's 
velocity vector to determine their effect on the line-of-sight 
angular rate (6) and the range rate (R). 

Figure 4-4 demonstrates the sign conventions used in this 
simulation. The trajectories which will be discussed later in 
Chapter Five use this same notation and sign convention. 

Figure 4-5 demonstrates in block diagram form the kine- 
matics. As noted in the figure, by-products of the kinematics 
are (1) the position of the missile and target, (2) the range and 
(3) the sign and magnitude of the range rate. 

By addition of a few fortran statements, adjusting the 
target's velocity vector, the target can be maneuvered as a 
function of the missile trajectory in range (evasion) or on a 
predetermined course (attack) . 

The fortran statements for the kinematics are listed below. 
C THIS SECTION DETERMINES THE EFFECT OF THE 
C PLANT OF THE KINEMATICS 
GAMDOT = X(1 ,1) 

GAMMA = GAMMA+GAMDOT*DT 
XMDOT = VM*COSF (GAMMA) 

YMDOT = VM*SIN,F (GAMMA) 

YRDOT = YTDOT-YMDOT 
XRDOT = YTDOT-YMDOT 

RDOT = ((YRDOT*SINF(SIGMA))+(XRDOT*COSF (SIGMA))) 
RDSIG = ((YRDOT*COSF(SIGMA))- (XRDOT*SINF (SIGMA))) 
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Figure 4-5. 



Block diagram of the Kinematics. 



RLOS = RLOS+RDOT*DT 
DSIG = RDSIG/RLOS 
SIGMA = SIGMA+DSIG*DT 
XMZ = XM Z+XM DOT *DT 
YMZ = YMZ+YMDOT*DT 
XTZ = XTZ+XTDOT*DT 
YTZ = YTZ+YTDOT*DT 
The Control 

In this simulation a control (A^) was selected through the 
use of program "OPCON" which minimizes the states X(k), ie, 
"bang-bang control" . Other forms of control are also avail- 
able such as placing a limit of the fuel or energy expended, 
placing a limit on the amount of acceleration allowed, or plac- 
ing a limit on the magnitude of the angle and angular rate that 
can be measured. These options are recommended for future 
study and would, of course, be necessary in making a more 
complete study of the missile systems. 

For this simulation the weighting vector (A 1 -) becomes a 
constant vector which is imposed on the difference between 
the best prediction of the states X*(k), and the deterministic 
input, DINP(k) to determine the control to be applied. 

Figure 4-6 is a discrete flow graph of the entire simulation. 
The fortran listing of the entire program may be found in Appendix 
One . 
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MISSILE DYNAMICS 







Figure 4-6. Discrete flow graph of missile simulation. 
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CHAPTER V 

SIMULATION RESULTS 

In making a study such as the evaluation of a missile guid- 
ance, it is quite difficult to decide which parameters to vary and 
which to hold constant. It is even more difficult to decide which 
curves hold the most meaning as a measure of success. Since 
the ultimate goal is to hit the target, trajectories and confirming 
print-out were used as a primary measure. Once a hit was ob- 
tained with a set of noise variances, the target trajectory was 
varied to see the effect on the missile. The noise was also 
varied to study the amount of noise that could be tolerated. 

The measurement noise was found to be the most critical 
value. A standard deviation of measurement noise of 0.1 radians/ 
sec was found to be the upper limit. Above this value the 
missile had some control, but the missile trajectories were far 
from desirable. Operation near the limit of measurement noise 
caused large ambiguities of the target information at the be- 
ginning of the flight. The target information improved as the 
range decreased, however, and in most cases the missile was 
able to capture the target. 

Curves of the normal acceleration of the missile, the filter 
gains, the range, and the control versus time were considered 
and are shown in the results where these values appear to be 
critical or of interest. 

This simulation assumed that the missile propulsion has 
already burned out and the missile is traveling at a constant 
velocity (V m ) . This simulation considers only the coplanar 
situation. The program could be easily expanded to two channels 
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Figure 5-1 INITIAL CONDITION OF TARGET 

Target velocity direction: A,B y C,D 
Missile: Same initial position for all runs 
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EVASION COURSES FOR TARGET 
Figure 5-2 TARGET TRAJECTORIES 
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Table 5-1 
INITIAL CONDITIONS 
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and crosscoupling could be considered. The trajectories may be 
considered as yaw or pitch maneuvers . If they are considered 
pitch maneuvers, the effects of altitude variations and gravity 
must be considered. Otherwise it is assumed that the missile 
has the same dynamic reaction to a change in either plane. The 
gravity consideration was not considered in this simulation. 

The target was made to under-go several maneuvers: (1) a 
turn, (2) a straight line course, and (3) a left/right/left turn 
(evasive course). Figure 5-1 is a graphic representation of the 
missile and target initial conditions where A, B, C, D represents 
the target velocity vector for the four cases. The numerical values 
of all the parameters involved are listed in Table 5-1 . The diff- 
erent maneuvers of the target are shown in Figure 5-2. A par- 
ticular computer run will be designated (A-l-2-0. 1-0.1). The 
"A" implies the initial conditions of the target (position and 
velocity direction). The first "l" implies that the magnitude 
of the target's velocity is 1000 ft/sec. The "2" implies that 
the target is flying a straight line course (maneuver two) . The 
last two numbers are the standard deviation of the excitation 
noise and the measurement noise respectively. Note that all 
angles in this program are given in radians . 

Figure 5-3 to Figure 5-2'6 are presented below to display 
the filtered missile capabilities. Listed on each figure are the 
range at which the program was terminated (R < 100 ft.) and the 
predicted miss distance. 
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MISSILE ACTIVE GUIDANCE SIMULATION 
A-l-2-0. 1-0.1 

Range (program termination) - 30. S ft. 
Predicted miss distance - 5.07 ft (right) 

Figure 5-5 
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MISSILE ACTIVE GUIDANCE SIMULATION 
A-l -2-0. 01-0. 01' 

Range (program termination) = 3.8 ft. 

predicted miss distance = 0.01 ft. (left) 

Figure 5-6 
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MISSILE ACTIVE GUIDANCE SIMULATION 
A-l-3-0 .1-0 *1 

Range (program termination) - 00. -St. 
Predicted miss distance - •f - *’* 

Figure 5-8 
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predicted miss distance = 38*9 ft. (left) 

Figure 5-10 




Cl 



ce 



os 

<s» 



vn 

Co 

c: 






8 iy 






kjCkj 









(Units = feet) 



K T Q Q 



~ n ^ 
. t — nu 



JNITS/IiSCH. 



! 1 if- 

iv^L 



B 



.1 ~ i ~a.l 



■id. 






GUIDANCE 



sinus 



or 

.ill 



T ^i 

1 i 1 



N 



53 



005 010 015 020 025 059 035 



Y 




MISSILE 






..HIT 



< 




target 



MISSILE ACTIVE GUIDANCE SIMULATION 
B-l~2-0.01-0.01 

Range (program termination) - 33.2 ft. 
Predicted miss distance - 0.05 ft. 



Figure 5—11 
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Predicted miss distance = 40.7 ft. (left) 

Figure 5-12 



809 095 810 915 818 815 850 

X 



* -SCALE " 5.805+93 UNITS/ INCH. 
Y- SCALE « 5.90E+83 UNITS/INCH. 



(units = feet) 



MISSILE AC'TIUE GUIDANCE SIMULATION 

R-1 -9-9.1 -0.1 55 



000 005 010 015 0-20 025 033 035 



Y 



MISSILE 




TARGET 



Range (program termination ) = 67.5 feet 
Predicted miss distance = 4.27 ft (right) 
Figure 5-13 
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Figure 5-14 
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Figure 5-16 
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APPENDIX I 



THROUGHOUT THIS THESIS COMPUTER PROGRAMS WERE USED. MADE 
REFERENCE TO OR EXPLAINED. THIS APPENDIX WILL LIST- THESE PRO- 
GRAMS AND GIVE A BRIEF EXPLANATION OF THE DATA REQUIRED. 

PROGRAM FILTER. 

THIS PROGRAM CALCULATES THE OPTIMUM STEADY-STATE GAINS 
G, FOR THE DIGITAL FILTER. THE DATA CARDS ARE EXPLAINED IN 
THE INITIAL COMMENT CARDS OF THE PROGRAM. 

PROGRAM FILTER 

C D1 ORDER OF SYSTEM IN 12 FORMAT 

C D2 SAMPLING INTERVAL IN 8F10.0 FORMAT 

C D3 F MATRIX BY ROWS IN 8F10.0 FORMAT 

C D4 D MATRIX BY COLUMN IN 8F10.0 FORMAT 

C D5 VARIANCE OF EXCITATION NOISE SIGWSQ IN 8F10.6 FORMAT 

C D6 R COVARIANCE MATRIX OF MEASUREMENT NOISE IN 8F10.6 FORMAT 

C PHI SYSTEM TRANSITION MATRIX 

C DEL DISTRIBUTION MATRIX 

C G OPTIMUM GAIN MATRIX 

C H OBSERVABLE MATRIX 

C P BEST ESTIMATE OF ERROR COVARIANCE MATRIX 

C Q EXCITATION NOISE COVARIANCE MATRIX 

DIMENSION P ( 12.12) »Q( 12*12 ) .H( 12*12) *R( 12*12 ) »G( 12*12) .PH IT ( 12. 12 
1 .PHI ( 12.12) *DEL ( 12 ) ,DELDELT( 12* 12) . DELS ( 1 2 » 1 2 ) . DELST ( 12.12 ) . 

2X( 1 2 ) . XHAT‘( 12 ) ♦ YHAT ( 12 ) .PNEW( 12.12 ) 

READ 33. N 
33 FORMAT (12) 

READ 2.DT 
2 FORMAT ( 8F10 • 0 ) 

PRINT 1003 
1003 FORMAT ( 1H 1 ) 

CALL PHIDEL(PHI .DEL.N.DT) 

DO 1001 LL= 1 .4 
READ 20. SIGWSQ 
20 FORMAT ( FI 0. 6 ) 

DO 1001 LK= 1 . 4 

C INIALIZE MATRICES FOR EACH RUN 

DO 10 I = 1 » N 
DO 10 J= 1 ♦ N 
PNEW ( I . J) =0.0 
P( I »J)=0.0 
G ( I ♦ J ) =0. 0 
Q( I .J) =0.0 
DELDELT ( I , J ) =0 .0 
RFLS H tr-h =0 *0 
10 DELST ( I . J ) = 0 . 0 
READ 2001 ,R( 1.1 ) 
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PROGRAM FILTER CONTINUED ! 

. i 

t 

2001 FORMAT ( 8F 10 #6 ) 

DO 30 1=1 .N 

30 DELS! I * 1 ) =DEL ( T ) 

CALL TRANS(DELS*N.N»DELST) 

CALL PROD ( DFLS» OFLST *N ,N*N,DELDFLT ) 

DO AO I =1 *N 
DO AO J=1 *N 

AO 0( I *J)=SIGWSQ*DELDELT( I * J ) \ 

SIGW=SQRTF(SIGWSQ) 

PRINT 100A * S I GW 

1 00 A FORMAT!/ 10X*5HSIGW= /.E20,8> 

P ( 1 * 1 ) = *02 l 

P ( 1 # 2 ) * 0*0 | 

P ( 2 * 1 ) = 0*0 | 

P!2*2>»1.0 t 

PRINT 2002* R ( 1 • 1 ) • ( ( 0! I ♦ J > . J = 1 *N ) . I =1 » N ) • ( ( P ( 1 ♦ J ) • J* 1 *N ) ♦ I = 1 ,N) 

2002 FORMAT! ///20X.8HR! 1*1)= • F10 .6/20X • 8H0{ I » J ) = »4F10* 6/20X *8HP ( 1 * J) » 

1 .AF10.6) i 

hi i * l ) = 1 • ;• 

H( 1 *2)30*0 

PRINT 700 f, 

700 FORMAT ( / /5X * 3HG1 1 # 7X ♦ 3HG12 *7X * 3HG2 1 *7X * 3HG22 *7X*3HP11»7X* 3HP12 * 7X * V 
1 3HP21 *7X*3HP22/> ( 

DO 1000 KK= 1 * AO | 

CALL GP ( H *PH I *P*Q*R*2*1 *G*PNEWJ 

PRINT 800! (G! I.J).J=1*N)*I=1.N) ♦( (PNEW! I ♦ J ) * J = 1 *N ) • I *1 * N ) 

800 FORMAT (10F10.5) 

DO 11 I = 1 *N i! 

DO 11 J=1 *N 

11 P ! I ♦ J ) =PNEW ( I *J> i 

1000 CONTINUE* f 

PRINT 1005 

1005 FORMAT! 1H1) J 

1001 CONTINUE f 

END f 

i 



SUBROUTINE GP ( H*PHI *P*Q*R.KN»KP*G*PNEW) 

DIMENSION H!12*12)*PHI!12.12).P(12*12)»0(12»12)*R(12»12).G! 12.12)* •' 
1 PNEW (12*12) »HT( 12*12) *TV1(12*12)*TV2(12*12) 

CALL TRANS!H»KP*KN»HTI i 

CALL PROD (P*HT*K.N»K.N»KP*TV1) 8 

CALL PROD ( H » TV1 *KP*KN*KP*TV2) 

CALL ADD(TV2*RtKP*KP*TVl ) 

CALL REC I P! KP t .0000000000001 *TVl *TV2 *KER ) 

I F ( KER-2 ) 101*110*101 1 

110 PRINT 111 S’ 

111 FORMAT (5HKER = 2) ’• 

101 CALL PR0D(HT*TV2»KN*KP*KP*TV1 ) j; 

CALL PROD (P*TV1*KN*<N*KP»G) 4 

gALL ) ) 

CALL PROD! G*TV1 *KN*KP*KN»TV2 ) j 
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PROGRAM FILTER CONTINUED 



DO 102 1=1, KN 
DO 102 J= 1 » KN 
102 TV2 ( I , J)=-TV2 ( I ,J) 

CALL ADD(P»TV2»KN»KN»TV1) 

CALL PRODtPHI ,TV1»KN,KN,KN,TV2) 
CALL TRANS ( PHI *KN,KN»TV1 ) 

CALL PR0D(TV2*TV1,KN,KN,KN,PNEW> 
CALL ADD(PNEW»Q*KN»KN»PNEW) 

END 

SUBROUTINE TRANS ( A ,N *M ,C ) 
DIMENSION A(12,12)»C(12*12) 

DO 153 I = 1 »N 
DO 153 J=1,M 
153 C ( J * I ) = A( I * J ) 

END 



SUBROUTINE R EC I P ( N , EP , A , X » KER ) 
DIMENSION A(12,12),X(12,12) 

DO 1 1=1, N 
DO 1 J= 1 , N 

1 X ( I > J)=0.0 
DO 2 K= 1 » N 

2 X ( K * K ) = 1 • 0 

10 DO 34 L=1,N 
KP = 0 

Z = 0.0 

DO 12 K=L,N 

IF(Z-ABSF(A 4 (K,L) ) )1 1,1 2, 12 

11 Z= ABSF ( A ( K » L ) ) 

KP = K 

12 CONTINUE 

IF (L-KP) 13,20,20 

13 DO 14 J =L , N 
Z=A ( L ♦ J ) 

A ( L » J ) = A ( KP » J ) 

14 A ( KP ♦ J ) =Z 
DO 15 J = 1 , N 
Z=X ( L , J ) 

X ( L , J ) =X ( KP ♦ J ) 

15 X ( KP , J ) =Z 

20 IF(ABSF( A(L,L) ) -EP ) 50 , 50 , 30 

30 I F ( L-N >31,34,34 

31 LP 1 =L+ 1 

DO 36 K=LP 1 , N 
I F ( A ( K , L ) ) 32,36,3 2 

32 RATIO=A(K,L)/A(L,L) 

DO 33 J=LP 1 » N 

33 A(K,J)=A(K»J)-RATIO*A(L,J) 

DO 35 

35 X(K,J)=X(K»J)-RATI 0*X ( L , J ) 

36 CONTINUE 
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34 CONTINUE 

40 DO 43 I =1 »N 
1 1 =N+1-I 

DO 43 J=1 »N 

S = 0.0 

I F ( I I-N)41,43.43 

41 IIP1=II+1 

DO 42 K= I I PI »N 

42 S=S+A ( I I ,K)*X( K.J) 

43 X(II»J)=(X(II»J)-S)/A(II,II) 
KER = 1 

RETURN 
50 KER = 2 
END 



SUBROUTINE PH I DEL ( PH I . DEL * N . DT ) 

C VALID ONLY FOR A CONSTANT F MATRIX 

DIMENSION F( 12*12) .PHI ( 12* 12 ) * TERM! 12* 12 ) *WORM( 12* 12 ) 

1. DEU 12) »DELM( 1 2 . 12 ) . TELM ( 12 . 12 ) . DELP( 12*121*0(12) 

READ1 * ( (F(IR*IC)*IC S 1*N)*IR=1*N) 

1 FORMAT ((8F10.0M 
READ 1 ( D( I ) * 1 = 1. N) 

1003 PRINT 399*DT.( (F(IR. IC).IC*1*N) *IR«1*N) 

399 FORMAT (///3HDT=*1F5.3///*7HF(I , J) */ , ( < 8F8.2 ) ) ) 

PRINT 3991 (D( I ) .1=1 *N) 

3991 FORMAT ( ////5HD( I ) =/( 8F8*2) ) 

NF I NAL = 1 
TM=0.0 

DO 400 I R= 1 » N 

DO 400 IC - 1 »N 

TERM(IR. 10=0.0 

WORM( IR.IC) =0.0 

TERM (IR»IR)=1.0 

TELM( IR.IC) =TERM( IR. IC)*DT 

DELP( IR.IC) =TELM( IR. IC) 

DELM( IR. 10=0.0 
DEL ( I R ) =0 .0 

400 PHI ( IR. IC)=TERM( IR.IC) 

4 TM=1.0+TM 

DO 500 IR=1.N 

DO 500 IC=1.N 

DO 500 JN=1.N 

DELM( IR. IC)=DELM( IR, IC)-TELM( IR. JN)*F( JN, IC) *DT/( TM+1.0 ) 
500 WORM( IR.IC) =TERM( I R» JN ) *F ( JN * I C ) *DT/TM+WORM( IR.IC ) 

DO 401 IR=1,N 
DO 401 IC= 1 »N 
TERM( IR. IC)=WORM( IR, IC) 

TELM( IR.IC) =DELM( I R. IC) 

DELP( IR.IC)=DELP( IR, IC)+TELM( IR.IC) 

r*M I ( leirpFM * ip» io^termi ir» io 

DELM( IR. IC) =0.0 

401 WORM (IR.IC)=0.0 
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ABC=0. 0 
DO 21 = 1, N 
DO 2 J= 1 »N 
AA=TERM( I ,J) 

AB=ABSF { AA) 

I F ( ABC-AB ) 3 , 3 , 2 
3 ABC=AB 
2 CONTINUE 

IF (0,00000000 5- ABC ) 5 ,5,6 

5 GO TO 4 

6 PRINT 502, ( (PHI ( I R , I C ) ,IC=1,N) ,IR=1,N) 
50 2 FORMAT ( ////9X,8HPH I ( I »J )/// (6E20.8 ) ) 

DO 600 1=1, N 
DO 600 K= 1 , N 
DO 600 J= 1 , N 

600 DEL ( I ) =DEL ( I )+PHI ( I,J)*DELP(J,K)*D(K) 
PRINT 503 < DEL ( I ) , I = 1 , N ) 

50 3 FORMAT ( // // 9X , 6HDEL ( I ) // ( 6E20 . 8 ) // ) 

END “ 



SUBROUTINE ADD (A,B»N»M,C) 

DIMENSION A(12»12)»B(12»12)»C(12»12) 
DO 152 1 = 1, N 
DO 152 J= 1 * M 

152 C ( I , J ) = A( I , J) + B ( I , J ) 

END 



SUBROUTINE PROD ( A , B , N , M, L , C ) 

DIMENSION A ( 1 2 » 1 2 ) ,B( 12,12 ) ,C( 12,12) * 

DO 151 1=1, N 

DO 151 J = 1 » L 
C(I,J) =0. 

DO 151 K = 1 , M 

151 C ( I , J ) = C(I,J) + A( I ,K)*B(K,J) 

END . 

END 



PROGRAM OPTCON 
DATA CARDS 

D1 CASE THE COMPUTER IS TO EXAMINE. 

CASE 1 Q = I, R = 0 , Q* = 0 
CASE 2Q=0, R=1 » Q* = 0 
CASE 3Q= I, R = 1* Q* EMPLOYED 
D2 N = ORDER OF SYSTEM, NSTAGE = NO. OF ITERATIONS 
D3 Q-WEIGHTING MATRIX 
R4 F MATRIX 
D5 D VECTOR 
D6 DT-SAMPLE PERIOD 
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PROGRAM OPTCON 

DIMENSION PH 1(12. 12). PS 1(12. 12). P(12. 12) .DEL ( 12) .AT ( 12) * 

1GM( 12.12) .0(12.12) »X( 900). IT I TLE( 1 2 ) . FM( 16 ) . EM(1 2 ) »HM( 12. 12 ) ♦ 

1 Y1 (900 ) .Y2 (900 ) . Y 3 ( 900) 

C THIS PROGRAM UTILIZES A COST FUNCTION* J ( N ) =MINIMUM( SUM X( N) T*Q*X( N) + 
C SUM R*U ( N- 1 ) **2 ) • AN UNLIMITED NUMBER OF ITERATIONS MAY BE MADE AT 
C A COMPUTATION RATE OF 2000 PER MINUTE AFTER THE PROGRAM HAS BEEN 
C COMPILED. THE OUTPUT OF THIS PROGRAM IS THE FEEDBACK GAIN MATRIX* 

C A TRANSPOSE. THE FOLLOWING RECURSIVE EQUATIONS WERE DERIVED USING 
C DYNAMIC PROGRAMMING. 

C AT( K)=-(DELT*P( K-l )*PHI ) / ( DELT*P ( K- 1 ) *DEL+R ) (1) 

C PSI (K) =PHI+DEL*AT(K) (2) 

C P(K)=PSIT(K)*P(K-1)*PSI (K)+Q+R*A<K)*AT(K) (3) 

C P(0 )=0,AT(0 )=0. PSI(0)=0 

C EQUATIONS 1* 2. AND 3 CONSTITUTE THIS PROGRAM 

C CALL IN DATA AND INITIALIZE 

DO 1111 1=1.3 
READ 30.KASE 
30 FORMAT ( ID 

READ 1 .N.NSTAGE 

1 FORMAT (8110) 

READ 2.R.DT 

2 FORMAT ( (4F20.0) ) 

READ 2 ( (Q( I .J) .J=1*N) .1=1* N) 

C PRINT THE DATA READ IN. 

PRINT 100 
100 FORMAT { 1 H 1 ) 

PRINT 32.KASE 
32 FORMAT (9X.5HKASE* .13) 

PRINT 3. N.NSTAGE 

3 FORMAT (//9X»2HN = *I3*20X* 7HNSTAGE* * 13) 

PRINT 4.R.DT 

4 FORMAT ( /9X*2HR=*F15. 11 *2X. 3HDT =. F15.ll) 

PRINT 9.( (Q( I »J).J=1.N).I=1*N) 

9 FORMAT ( /9X.7HQ( I,J)=/(2F15.11)) 

CALL PHIDEL(PHI.DEL.N.DT) 

DO 5 I =.l ♦ N 
DO 5 J= 1 *N 
GM ( I .J)=0.0 
EM ( I ) =0 .0 
FM ( I )=0.0 

5 P( I . J ) =Q( I .J) 

C INITIALIZE P(0) FOR CASE TWO 

I F ( KASE-2 ) 35.36*35 
36 P(l.l)=1.0 
P( 2.2 ) =1.0 
P(3,3)=1.0 
35 CONTINUE 
PRINT 19 

19 FORMAT ( //2X » 6HNSTAGE • 4X *7HAT ( 1 * 1 ) *3X *7HAT ( 1 * 2 ) *4X *6HP( 1*1)* 
1*X*6HP(1»2) * 4X * 6HP (2.1 ) .4X *6HP ( 2 *2 ) » 

C CALCULATE AT(J) 

DO 22 KK=1 .NSTAGE 
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DEN=0. 0 
D06 1=1, N 
D06 J= 1 ,N 

6 EM( I )=EM( I ) +DEL I J ) *P I J » I ) 

D08 1=1, N 

D07 J= 1 , N 

7 FMI I )=FM( I )+EM( J)*PHI I J,I ) 

8 DEN=DEN+EM ( I )*DEL( I ) 

DEN=-DEN-R 

DO 101 = 1, N 
AT { I ) = FM ( I ) /DEN 
FM { I ) = 0 .0 
10 EMI I ) =0 .0 

C CALCULATE PSI(K,I,J) 

DO 131 = 1, N 
DO 1 3 J= 1 »N 

13 PS I I I , J ) = PHI ( I , JJ+DELI I )*AT ( J ) 

C CALCULATE P(K,I,J) 

DO 15 I = 1 , N 
DO 15 J = 1 ,N 
DO 15 L = 1 , N 

15 GM ( I ,J)=GM( I ,J)+PSI I L » I ) *P I L » J ) 

DO 1 7 1 = 1 , N 
DO 17J= 1 ,N 
DO 16 L = 1 ♦ N 

16 HM( I , J ) =HM ( I , J ) +GM ( I ,L)*PSI (L,J) 

CASE 1 TERMINAL CONTROL, OMIT Q ( I * J ) IN EQUATION FOR PI I , J ) 
CASE2 MINIMUM FUEL OMIT Q ( I ♦ J ) IN EQUATION FOR PI I , J ) 

I F ( KASE-2 > 31,31,33 
31 P ( I ♦ J ) =HM ( I , J ) +R*AT I I ) *AT I J ) 

GO TO 17 

C CASE THREE MINIMIZATION OF TIME AND FUEL 

3 3 P ( I ♦ J ) =HM ( I , J ) +Q ( I , J ) +R*AT ( I ) * A T ( J ) 

17 HM( I ,J ) =0.0 
DO 18 I.= 1,N 
DO 1 8 J= 1 ♦ N 

18 GM( I , J ) =0.0 

PRINT 20,KK,AT(1),ATI2),P(1,1),P(1,2),P(2,1),P<2,2) 

20 FORMAT ( 15, { 6F10.6) ) 

22 CONTINUE 
1111 CONTINUE 
END 

SUBROUTINE PH I DEL ( PH I ♦ DEL , N , DT ) 

C VALID ONLY FOR A CONSTANT F MATRIX 

DIMENSION FI12, 12), PH 1(12,12) , TERM ( 12,12 ) , WORM! 12,12 ) 

1 , DELI 12 ) »DELM I 12*12)»TELM(12,12) ,DELPI 12,12), D(12) 

1 FORMAT (I8F10.0)) 

READ1 , I (FI I R , I C ) , I C= 1 » N ) ,IR = 1,N) 

READ 1 (DU) , 1 = 1 ,N ) 

1003 PRINT 399, DT , I I FI I R, IC ) , IC = 1 ,N) , I R = 1 ,N) 

399 FORMAT I ///3HDT=,1F5.3/ , 7HF I I , J ) =/ , ( I 2 F 8 . 2 ) ) ) 
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PRINT 3991 (DU) , I =1 *N> 

3991 FORMAT ( / 5HD( !)=»/( 2F8.2 ) ) 

NF I NAL= 1 
TM=0.0 

DO 400 IR=1 »N 

DO 400 IC=1 »N 

TERMUR.IC) =0.0 

WORMUR, IC)=0.0 

TERM ( IR»IR)=1.0 

TELM ( IR.IC) =TERM( IR. IC)*DT 

DELP( IR.IC)=TELM( IR. IC) 

DELM ( I R * I C ) =0 .0 
DEL ( I R ) =0 • 0 

400 PHI ( IR.IC)=TERMt IR.IC) 

4 TM= 1 • O+TM 

DO 500 IR=1 »N 

DO 500 I C=1 »N 

DO 500 JN=1.N 

DELM( IR.IC) =DELM( IR. I C )-TELM( I R » JN) *F( JN. IC ) *DT/ ( TM+1.0 ) 
500 WORMUR, IC)=TERM( I R , JN ) *F ( JN. I C ) *DT/TM+WORM< IR.IC) 

DO 401 I R = 1 » N 
DO 401 IC = 1 ,N 
TERM! IR.IC) =WORM( IR.IC) 

TELM( IR.IC) =DELMC I R. IC) 

DELPI IR.IC >=DELP( I R, IC)+TELM( IR.IC) 

PHI ( IR, IC)=PHI ( IR, IC)+TERM( IR, IC) 

DELM( IR.IC) =0.0 

401 WORM! IR. 10=0.0 
ABC=0. 0 

DO 21 = 1. N 
DO 2 J= 1 ,N* 

AA = TERM( I ,J) 

AB=AB5F(AA) 

I F ( ABC-AB ) 3 .3 ,2 
3 ABC =AB 
2 CONTINUE 

IF(0.00 0000005 -ABC)5*5»6 

5 GO TO 4 

6 PRINT 11* ( (PHI ( I »J)»J=1»N) ,I=1,N) 

11 FORMAT (/9X.9HPH I ( I ,J )=/( 2F 15. 1 1 ) ) 

DO 600 1=1. N 

DO 600 K= 1 . N 
DO 600 J = 1 * N 

600 DEL ( I ) =DFL ( I )+PHI ( I . J ) *DELP ( J ,K ) *D ( K ) 

PRINT 12 » CDEL( I ) ♦ I =1 ,N ) 

12 FORMAT (9X,7HDEL( I ) =/ ( 1F15. 11 ) ) 

END 

END 
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PROGRAM MISSILE 

C THIS PROGRAM SIMULATES THE MISSILE CONTROL , DYNAM I CS * AND KINEMATE 

THE KALMAN FILTER AND OPTIMAL CONTROL IS APPLIED TO THE CONTROL 
C OF THE PLANT. 

C THIS PROGRAM USES VARIABLE FILTER GAINS 

C 

C THIS PROGRAM CLOSES ThE LOOP ON THE OPTIMUM FILTER-CONTROLLER 

C PROBLEMS. IT ASSUMES THAT STABLE CONTROLLER GAIN MATRIX HAS 

C BEEN COMPUTED ON THE BASIS OF DESIRED RESPONSE AND CALCULATES 

C THE FILTER GAIN MATRIX EACH SAMPLE ON THE BASIS 

C OF THE STATISTICAL PROPERTIES OF THE ANTICIPATED RANDOM DISTURBA! 

C 

C THE PROGRAM SOLVES THE FOLLOWING EQUATIONS 

C Y(K)=H*X(K) 

C Z ( K ) = Y ( K ) +V ( K ) 

C XS(K) = ( I-GH)PHI*XSI.<-1 ) + ( I-GH)*DEL*AT(XS(K-1 )-DINP( K-l) )+G*Z(K)‘ 

C X(K+1)=PH!*X(;<) +DEL*Vv ( ,< ) tDEL*AT* ( XS { K ) -D I NP ( K ) ) 

C WHERE V IS MEASUREMENT NOISE*',*.' IS THE RANDOM D I STUR8ANCE » AND 

C D I NP IS THE DETERMINISTIC INPJT 

C 



DATA CARDS. 

D 1 N = ORDER. OF SYSTEM* DT = SAMPLE PERIOD 
D2 GAMMA = FLIGHT PATH ANGLE 
SIGMA = LOS ANGLE 

GAMDOT = FLIGHT PATH ANGULAR -RATE 
DSIG = LCS ANGULAR RATE 

D3 RLOS = RANGE 

RZ = REQUIRED TARGET RANGE FOR HIT 
XTZ = INITIAL X-COORDINATE OF TARGET 

XTDOT = INITIAL VELOCITY IN X DIRECTION FOR THE TARGET 
YTZ = INITIAL Y-COORDINATE OF TARGET 

YTDOT = INITIAL VELOCITY IN Y DIRECTION FOR THE TARGET 
DA INITIAL CONDITIONS OF THE STATE VECTOR 
D5 XMZ INITIAL X-COORDINATE OF MISSILE 
YMZ INITIAL Y-COORDINATE OF MISSILE 
VM MISSILE SPEED 
D6 P-BEST ESTIMATE OF ERROR - 
D7 F MATRIX 
D8 D VECTOR 

D9 GRAPH TITLE lINE ONE 

DiO AT FEEDBACK GAIN MATRIX FOR CONTROLLER 
DU GRAPH TITLE LINE TWO 

D 1 2 S I GW = STANDARD DEVIATION OF EXCITATION NOISE 
SIGV = STANDARD DEVIATION OF MEASUREMENT NOISE 

DIMENSION X( 12*12) *XS( 12*12) »SIGV( 12) »Y ( 12*12 ) *Z( 12* 12) *PHI ( 12*12 

1 » DEL ( ’2 >12) *H( 12*12) * A T (12*12) » G ( 1 2 » 1 2 ) , TEMPI ( 12,12) ,TEMP2( 12,12) 

2 TEMP 3 ( 12*12) * TEMP A ( 1 2 * 1 2 ) *TELP( 12*12), TELP 1(12, 12) >TELP2( 12*12) * 

3 V ( 12,12) ,DINP( 12*12) * IT (12) »TME(900) ,YT(900) , XT (900) 

A » X M ( 9 0 0 ) , Y M ( 9 0 0 ) ,P(12*12)*R(l2»12) 
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5 * DELTR ( 12,12 5 » DELS ( 1 2 5 , Q < 1 2 , 1 2 5 
6 , SDOT ( 900 5 , D I FF ( 90 0 ! 

C INPUT SOME CONSTANTS, AS NOTED 

READ 2 » N > DT 

READ 1 , GAMMA , S I GM.A » GAMDGT » DS I G 
READ 1 ,RLOS,R2,XTZ , XT DOT ,YTZ ,YTDOT 
READ 1 , CXC I .1 ) * 1 = 1 »N> 

READ 1 » XMZ , YMZ » VM 
3 FORMAT (6AS) 

2 FORMAT C 15 , IF 10. 0) 

1 FORMAT (( 8F10.0 ) 5 

READ 1 , { C PC I , J) »j=l »N 5 , 1 = 1 ,N 5 
CALL PHIDELCPHI , DELS , N , DT 5 
READ 3, ( I T t I 5 » 1=1,6) 

READ 1 » ( AT ( 1 » I ) » 1 = 1 ».\ i 
READ 3 » ( I T ( I 5 » I =7 , 12 5 
READ 1 ,SIGW»SIGV< 1 5 
PRINT 144 

144 FORMAT ( 4X , 8HXM I SS I LE , o X , SriXTARGET , SX » SriYM ISS I LE » 8X » 8HYT ARGE T ) 
DO 6000 I = 1 , N 
6000 DFLC I , 1 5=DFLS( I 5 

CALL TRANS ( DEL , N »N , DELTR ! 

CALL PROD ( DEL , DELTR > N / \ > N » 0 5 
DO 2 00 j=;,,\ 

DO 200 1=1, N 

200 q ( i , j ) =c ( : , j ) *s i g.--*s : gw 
R( i , i) = s:gv( : 5*s:gvc i ; 

P(1,1)=r;:->15 

C AN IS SYSTEM ORDER, .<? IS Tr.E NUMBER OF OBSERVABLES. 

KN =N 
KP = 1 

H { 1 , 1 ) = 1 . 0 

C SETTING DINP=0.C AT TrilS POINT INSURES THAT NO DETERMINISTIC 

C INPUTS EXIST PRIOR TO TIME = ZERO 

DO 131 I = 1 » ,< N 
131 XSC I, 15=0.0 

NUN IF =1220703125 
C FILTER INITIAL CONDITIONS 

CAwL PROD ! H , X » ;<? , AN •. :<? > Y 5 
DO 12 !=!,;<? 

CALL RN DEV ! NUN IF, DEV ) 

12 VC I ,15 =SIGV< I >*DEV 
CALL ADD! Y,V,AP, 1*Z 5 

xs c i ,1 ; =zc i ,i ) 
k;<= i 
T = 0 ,0 
KK 1 = 0 

C THIS SECTION CALCULATES G, THE GAIN MATRIX 

/■* 

6009 CALL F I w 7 ER ( AN , AP , P.H I » 0 » R » H » P » G 5 
C 

C THIS SECTION CA i _C0_ATES XS, THE BEST ESTIMATE 



81 



n n n n o n 



PROGRAM MISSILE CONTINUED 



C OF THE STATE VECTOR 

CALL PRODt H,X »KP,XN, 1 »Y ) 

DO 10 1=1, KP 
CALL RNDEV(NUNIF,DEV) 

10 V ( I , 1 ) = S I GV ( I ) *pEV 
CALL ADD(Y*V»KP»1*Z) 

CALL PROD ( PH I » XS > KN » <N , 1 * T EM? 1 ) 

CALL PROD(G*H»KN»KP*.<NtTEMP2 J 
DO 11 1 = 1 5 KN 
DO 11 J = 1 » K N 

11 TEMP2( I,J)=-TFMP2( I, J) 

CALL PRODt TEMP 2 , TE MR 1 v < N , K N , 1 , TEMP 3 ) 

CALL ADD (TEMP 1 , TEMP 3 , KN » 1 , TEMP 3 ) 

CALL ADDt XS » D I NP , KN , 1 ,TELP ) 

CALL PRODt AT »TELP* 1*<N» 1, TEMPI ) 

CALL PRODt DEL * TEMPI »<N » 1 * 1', TEL? ) 

CALL PROD(TEMP2»TELP*KN*KN»l ,TELP1 ) 

CALL ADD( TELP »TELP1*KN» 1 .TELPl ) 

CALL PROD(G»Z»KN*<P* 1 » TELP2 ) 

CALL ADD ( TEMPS » TELP1 » .<N ♦ i » XS ) 

CALL A DDT XS,TELP2»KN » 1 »XS ) 

D I NP ( 1 , 1 ) =-DS I G 
D IMP (2,1) =-DDSIG 
SDOT ( KK ) =DS I G 
TME ( KK ) =T 

THIS SECTION UPDATES THE STATE VECTOR X 

CALL RNDEVt NUM F » DEV ) 

W=SIGV/*DEV 

CALL PROD { PH I > X » KN , KN , 1 , TEMP 1 ) 

DO 803 1=1 , KN 
803 TEMP2 (1,1) =W*DEl (1*1) 

CALL ADD ( XS , D I \'P »KN , 1 *TELP ) 

CALL PRODt AT ,TEl=* 1 » KN » 1 , TELPl ) 

CALL PROD (DEL»TEuPl»KN»l»l»TELP2) 

CALL ADD ( TEMP 1 , JE: j .?2 ,KN,i ,X ) 

CALL A D D ( X , T E L P 2 , N > 1 » X ) 

THIS SECTION DETERMINES THE EFFECT OF THE PLANT 
ON THE KINEMATICS 
DKS I G=DS I G 

DI FF ( KK ) =X ( 1 , 1 i -XS t 1,1) 

GAMDOT = X ( 1,1) 

GDOT(KK) =GAMDOT 
GAMMA =GAM M-.-rG,A' / DO T^DT 
XMDOT =VM*COSE ( GAMMA ) 

133 YMDOT = VM-"'S I ,\F ( GAMMA ) 

1 3 A YRDOT =YTDOT-YMDOT 
XRDOT =XTDOT-XMDOT 

ROOT = ! (YRDOT*S INF (SIGMA) ) + ( XRDOT*COS F ( S I GMA ) ) ) 
RDS I G= ( ( YRDOT *COSF( SIGMA) ) - £ XRDOT*S I NF ( S I GMA ) ) ) 

sa 
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RLOS =RLOS+RDOT*D7 
DS I G=RDS I C- 'RLOS 
SIGMA =SIoMAtDSIG*DT 
DOS I 6 = { DSIG-DKSIGi /D7 
XMZ = XMZ+XMDOT*OT 
YMZ = YMZ+YMDOT*DT 
XM t K !< ) -XMZ 
YM ( KK 5 = YMZ 
XTZ -XTZ+XTDOT*DT 
YTZ =YTZ+YTDGT*DT 
YT ( KK ) =YTZ 
X7 cco =xtz 
RL ( KK ) =RLOS 
PRINT 140 5 KK 

PRINT 141 »XM( KK ) jXT ( KK ) » YM i KK ) > Y T ( KK } 

140 FORMAT; 2X» 15) 

141 FORMAT ( 2X *4E16„7) 

KK=KKt 1 

7 = T + D T 

PRINT 190s DSIGsSIGMA 

190 FORMAT ( 6n DS IG= v IE 16 .7 > 2X *7h S I GMA = > 1 E 1 6 . 7 ) 
IF ( ROOT 5 146} 145 , 145 
149 PRINT 1 35 » ROOT 
GOT 0 6013 

146 I F { K K— 9 0 0 ) 6 C 1 2 v 6 0 1 2 > 6 0 1 3 

6012 IF (RLOS- ICO. ; 60 1 3 > 6C 1 3 , 60 1 1 
6011 GOTO 6009 

6013 PRINT _ 36 » RLOS 

1 3 FOR I' 1 * 6u kDG i = j.;1o. i 1 

136 FORMAT ( 6H RLOS- , IE 16.7 1 
A — A . A\r ( R D S I G / ^00 1 i 

tm i ss= { :rlos*ros:g‘, /root j*cosf{ aj 

PRINT 137 > 7 M I S S 

137 FOR. -'AT ( 7i-i TM I SS = •, 1 E 1 6 „ 7 ) 

KK-K.n-. 

MC -1 



0RA1.'( KKvXM 9 YM»MCvO»wA> I i •» 0 9 0 » 0 » 0 » 0 » 0 > 7 » S 1 



MC - 3 
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PROGRAM MISSILE CONTINUED 



DO 152 1=1, N 
152 C ( I * J ) = A ( I , J)+B( I , J) 
END 



SUBROUTINE F I LT ER ( N , KP ♦ PH I ,Q ♦ R ,H ,P ,G ) 

PHI SYSTEM TRANSITION MATRIX 
DEL DISTRIBUTION MATRIX 
G OPTIMUM GAIN MATRIX 
H OBSERVABLE MATRIX 

P BEST ESTIMATE OF ERROR COVARIANCE MATRIX 
Q EXCITATION NOISE COVARIANCE MATRIX 
DIMENSION P( 12, 12) »Q( 12,12 )•» H ( 12,12) »R( 12,12 ) »G( 12, 12) »PHIT ( 12»12 
1 , PHI ( 12, 12) , DEL ( 12 ) , DELDELT (12*12) ,DELS( 12*12) , DELST ( 12,12 ) » 
2PNEWI 12,12) 

CALL GP(H,PHI » P » 0, R » N » KP » G » PNEW ) 

DO 11 1=1, N 
DO 11 J=1,N 

11 P ( I *J) =PNEW( I , J) ' 

END 

SUBROUTINE GP ( H , PH I , P , 0 , R ♦ KN , KP , G , PNEW > 

DIMENSION H( 12*12 > *PHI ( 12*12) *P( 12*12) *Q( 12*12) *R( 12*121 *G( 12*12) 
1 PNEW ( 12*12), HT( 12, 12) *TV1(12*12)*TV2 (12*12) 

CALL TRANS(H»KP»KN»HT) 

CALL PROD (P»HT »KN»KN»KP,TV1 ) 

CALL PR0D(H,TV1 , KP , KN » KP , T V2 ) 

CALL ADD( TV2 ,R*KP,KP»TV1 ) 

CALL REC I P ( KP, .0000000000001 ,TVl»TV2»KER) 

I F ( KER-2 ) 101,110,101 I 

110 PRINT 111 

111 FORMAT ( 5HKER=2 ) 

101 CALL PROD (HT »TV2,KN,KP,KP»TV1) 

CALL PR0D(P,TV1,KN,KN,KP,G) 

CALL P ROD ( H , P , KP , KN , KN ♦ T V 1 ) 

CALL PROD ( G » T V 1 , KN , KP , KN » T V2 ) 

DO 102 1=1, KN 

DO 102 J= 1 * KN 

102 TV2 ( I , J ) =-T V2 ( I , J ) 

CALL ADD( P,TV2,KN,KN,TV1 ) 

CALL PRODtPHI , T Vl , KN , KN , KN , T V2 > 

CALL TRANS ( PH I ,KN,KN,TVl) 

CALL PROD( TV2,TV1,KN,KN,KN,PNEW) 

CALL ADD(PNEW,Q,KN,KN,PNEW) 

END 

SUBROUTINE TRANS ( A ,N ,M ,C ) 

DIMENSION A< 12, 12) ,C( 12 ,12 ) 

DO 153 I = 1 , N 
DO 153 J= 1 , M 
183 H = AH 

END 
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PROGRAM MISSILE CONTINUED 



SUBROUTINE TRANS ( A ,N ,M , C ) 
DIMENSION A(12*12),C(12.12) 

DO 153 I = 1 * N 
DO 153 J = 1 » M 
153 C C J * I ) = A{ I * J ) 

END 

SUBROUTINE R EC I P < N » EP , A , X . KER ) 
DIMENSION AU2.12) ,X(12.12) 

DO 1 1=1 .N 
DO 1 J*1,N 

1 X ( I , J) =0.0 
DO 2 K = 1 » N 

2 X(K.K)=1.0 

10 DO 34 L= 1 *N 
KP = 0 

2 = 0.0 

DO 12 K = L,N 

IF(Z-A6SF(A(K,L> )) 11.12.12 

1 1 Z=ABSF( A( K,L) ) 

KP = K 

12 CONTINUE 
IFIL-KP) 13,20.20 

13 DO 14 J=L,N 
Z = A ( L » J ) 

A ( L » J ) = A ( KP , J ) 

14 A ( KP , J ) =Z 
DO 15 J = 1 , N 
Z=X ( L. J ) 

X(L.J) =X(KP, J) 

15 X ( KP , J ) =Z 

20 IF(ABSF(A(L,L) >-EP >50.50.30 

30 IF(L-N) 31 ,34,34 

31 LP 1 =L+ 1 

DO 36 K = LP 1 , N 
IF(AIK.L) >32,36,32 

32 RATIO=A(K.L)/A(L»L) 

DO 33 J=LP1.N 

33 A(K.J) =A(K, J)-RATIO^A(L.J) 

DO 35 J=1 , N 

35 X(K,J)=X(K,J)-RATIO*X(L»J) 

36 CONTINUE 

34 CONTINUE 

40 DO 43 I = 1 . N 
I I =N+1-I 

DO 43 J=1*N 
S=0.0 

I F ( I I -N >41,43,43 

41 I I P 1 = I I + 1 

DO 42 K=I IP1 ,N 

42 s=s*a< u < K**n 

43 XI I I . J > =(X( I I , J)-S) /A( I I • I I ) 
KER = 1 

SS 
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PROGRAM MISSILE CONTINUED 



RETURN 
50 KER=2 
END 

SUBROUTINE PHIDEL ( PH I , DEL , N , D T > 

DIMENSION F(12»12)* PH 1(12*12) , TFRM ( 12*12) ,WORM( 12,12) 

1* DEL ( 12) ,DELM( 12,12) *TELM(12,]2) ,DELP( 12*12), D ( 12 ) 

VALID ONLY FOR A CONSTANT F MATRIX 
DT = SAMPLING INTERVAL 
NF I NAL= 1 

READ 1( (F( IR»IC),IC=1,N)»IR=1,N) 

READ 1 (D( I ) , 1 = 1 ,N) 

1 FORMAT ((8F10.0)) 

TM=0 • 0 

1003 PRINT 399*DT ,( ( F ( IR» IC) » I C = 1 » N ) » I R = 1 * N ) 

399 FORMAT ( / / /4H DT = * 1F5.3////,8H F ( I , J ) =/ , 4 ( 4E 1 6 . 7/ ) ) 

PRINT 3991 (D(I)*I=1»N) 

3991 FORMAT ( / / /6H D ( I ) = / * 4 ( 1 FI 6 . 7/ ) ) 

DO 400 I R = 1 , N 

DO 400 I C = 1 » N 

TERM ( I R * I C ) =0.0 

WORM ( I R » I C ) =0.0 

TERM ( I R, IR) =1 .0 

T E L M ( I R, IC ) =T ERM ( I R, IC )*DT 

DELP ( I R , IC ) =TELM( IR, IC ) 

D E L M ( IR , IC ) =0.0 
DEL ( IR ) =0.0 

400 PHI ( IR, IC)=TERM( IR, IC) 

4 TM= 1 . O+TM 

DO 500 I R = 1 , N 

DO 500 I C = 1 » N 

DO 500 JN = 1 , N 

DELM ( IR,IC)=DELM( IR, IC)-TELM( I R , JN ) *F < JN , I C ) *DT / ( T M+ 1 . 0 ) 
500 WORM ( IR, IC) =TERM( IR, JN)*F ( JN, I C ) *DT / TM+WORM ( IR, IC) 

DO 401 I R= 1 , N 
DO 401 I C = 1 , N 
TERM ( IR, IC) =WORM ( IR, IC ) 

TELM ( IR,IC)=DELM( IR, IC) 

DELP ( IR, IC) =DELP ( I R , I C ) +T ELM ( I R , I C ) 

PHI ( IR , IC)=PHI ( IR, I C ) +TERM ( IR , IC) 

DFLM ( IR,IC) =0.0 

401 WORMdR, 10=0.0 
ABC=0 • 0 

DO 2 1 = 1 ,N 
DO 2 J= 1 ,N 
AA = TERM ( I , J ) 

AB=ABSF ( AA) 

I F ( ABC-AB ) 3 * 3 , 2 
3 ABC=A8 

2 CONTINUE 

I F ( 0.000000005-ABC ) 5 ,5 ,6 

H (EjE) (C) 4 

6 PRINT 502 , ( < PHI ( IR , IC) » IC=1,N) , IR=1 ,N) 

502 FORMAT ( //9X,8HPHI ( I,J)///((4El6.7))) 
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PROGRAM MISSILE CONTINUED 



88 888 i:l:H 

DO 600 J= 1 » N 

600 DFL ( I ) =DEL ( I J+PHI ( I , J ) *DELP ( J . K ) *D ( K> 
PRINT 503 ( DEL ( I ) • I = 1 * N ) 

50 3 FORMAT ( // /9X * 6HDEL ( I)//(8F15.Q>//) 
RETURN 
END 
END 



PROGRAM MISSILE CONTINUED 



C EVASION 

YTDOT =YTD0T+10» , 

C EVASION LEFT-RIGHT-LEFT 

IF(<K-3D)320»321,321 

320 YT DOT = YTDOT +20# 

GOTO 323 

321 YT DOT = YTDOT -20* , 

323 CONTINUE 

I F ( YTZ-5000 . ) 324,324,325 

324 YT DOT = YTDOT +30. 

IF (YT2-2000. ) 326,326,325 
326 YTDOT = 0*0 

325 CONTINUE 

C END OF EVASION 
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APPENDIX II 



The missile guidance simulation is a transfer function relating 
( y ) and ( o ) . Since the objective of this study is to examine the 
effect of a digital filter, an arbitrary transfer function was se- 
lected which could represent any present day missile. This 
transfer function contains a time delay due to the radar, a fil- 
tering action due to the control (hydraulics etc.), and a second 
order system representing the missile dynamics. 

o/y = A . B . C (A-l) 

S + T r S + T c S 2 + DS + E 

As mentioned above the intention here is to demonstrate the 
capability of the digital filter in the presence of a large mag- 
nitude of noise. The values of the gains and time constants were 
selected for equation (A-l) after examining several different 
missile systems which use proportional navigation. Refer to 
[13] for further elaboration on determining the F and D matrices 
for a particular system. 
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